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ABSTRACT 

A useful crude approximation for Abelian functions is developed and applied to orbits. 
The bound orbits in the power-law potentials Ar~ a take the simple form (£/r) k — 
1 + ecos(m0), where k = 2 — a > and £ and e are generalisations of the semi- 
latus-rectum and the eccentricity, m is given as a function of 'eccentricity'. For nearly 
circular orbits m is y/k, while the above orbit becomes exact at the energy of escape 
where e is one and m is k. Orbits in the logarithmic potential that gives rise to 
a constant circular velocity are derived via the limit a — » 0. For such orbits, r 2 
vibrates almost harmonically whatever the 'eccentricity'. Unbound orbits in power- 
law potentials are given in an appendix. The transformation of orbits in one potential 
to give orbits in a different potential is used to determine orbits in potentials that 
are positive powers of r. These transformations are extended to form a group which 
associates orbits in sets of six potentials, e.g. there are corresponding orbits in the 
potentials proportional to r, r~ 3 , r -3 , r~ 6 , r~ 3 and r 4 . A degeneracy reduces this to 
three, which are r -1 , r 2 , and r -4 for the Keplerian case. A generalisation of this group 
includes the isochrone with the Kepler set. 
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1 INTRODUCTION 

Since schooldays when we encountered the rigid pendu- 
lum, most of us have been frustrated by our inability to 
integrate in elementary terms Abelian expressions of the 
form J [S(u)]~z du, where S(u) has simple zeros at u a and 
Up ^ u a but is not quadratic. In practice S usually depends 
linearly or quadratically on parameters which we shall call 
e and h, and its zeros u p and u a depend on e and h often in 
quite complicated ways. In Appendix A we relieve this frus- 
tration by showing how for each pair of u a and u p , S may 
be replaced at lowest order by a quadratic function with 
the same zeros and the integral may be approximately 
evaluated parametrically via perturbation theory. 

We do not have to solve S(e, h, u) for its zeros u p and u a . 
Instead we regard u p and u a as parameters and then easily 
find the s(u a ,u p ) and h(u a ,u p ) to which they correspond. 
The process of replacing 5 by a different quadratic function 
for each pair of zeros u p and u a we call quadrating (after 
the old verb 'to quadrate' which means 'to make square'). 
Surely using 'quadrate' to mean 'to make quadratic' is not 
too great an extension! As the simple, though crude, method 
developed can be applied to a far wider class of problems 
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than those encountered here, we have mentioned it first in 
the introduction. 

Orbits of the general form 



(i/r) 



1 + e cos (mcf 



(1) 



have a long history. Newton in Principia l|l687t) showed that 
orbits of this form with k = 1 occurred when the central 
force was an inverse square law supplemented by an in- 
verse cube force. His celebrated theorem on revolving orbits 
demonstrates that if r = r((j>) is an orbit of angular mo- 
mentum h under any central force F(r)r, then r = r(m<j>) is 
an orbit of angular momentum mh under the central force 
[F — (rri 2 — l) ft 2 r~ J ] f . Newton also pointed out that r(t) 
was the same for both orbits, so, if the new orbit were viewed 
from axes revolving at the rate (m — l)<j>, then the two or- 
bits would have the same shape. But notice that (m — 1)0 
is not a constant rotation rate, but speeds up when r is 
small and slows down when r is large. When we view the 
orbit from axes that rotate uniformly at the same mean 
rate (m — 1) < (j> >j the orbits can have very different 
shapes involving figures of eight fo r the more eccentric ones 

jLvnden-Bell fe Lvnden-Bel]||l995l'). 

Recently in a fine paper, IStruckl (|2006l ) showed that or- 
bits of moderate or low eccentricity in logarithmic or power- 
law potentials with or without cores were well approximated 
by analytic orbits of the form ([l]). His approximate orbits 
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are surprisingly accurate. Struc k was, in part, stimul a ted to 
find this result by a paper by iTouma fc Tremainel (|f997l ) 
that demonstrated the richness of the resonance s in th e per- 
turbation theory of these systems. IValluri et all (|2005l ) have 
emphasised that the apsidal precession found in non-inverse- 
square orbits is significantly dependent on the eccentricity 
of the orbit involved. 

Surprisingly, we have been led to orbits of the form fT} 
by looking at orbits of extreme eccentricity e = 1, where 
Struck's metho ds did not give accura t e res u lts. Standard 
works on orbits. | Boccaletti fc Pucaccol l|l996t ) , IContopoulosI 
l|2002l) and lBinnev fe Tremainel (Il987l ). do not point out that 
the orbits of zero energy in power-law potentials can be ex- 
actly solved analytically. The same variables can be used to 
solve the nearly circular orbits. Since both highly eccentric 
and small eccentricity orbits can be so solved, it would be 
surprising if there were not a good approximation, based 
on the same variables, that interpolated between e = 1 and 
e = 0. Struck's methods do this well for small and moderate 
eccentricities. Here we show that all orbits are well approx- 
imated by analytic orbits of the form {TJ for ^ e ^ 1. 

Here £ and e are generalisations of the semi-latus- 
rectum and the eccentricity. The potential is Ar~ a = Ar k ~ 2 \ 
this parameterisation, using fc rather than a, is chosen in this 
section to simplify equations ((2} and (|3} below. If r a and r p 
are the apocentric and pericentric distances the generalised 
eccentricity is given by 



(2) 



+ r a 



and the generalised semi-latus-rectum is given by £ where 



5" 



p f a 



(3) 



We note that for the Kepler case fc — 1 and the above 
formulae all reduce to the usual ones. We write £ — Lr c , 
where r c (h) is the radius of the circular orbit of angular mo- 
mentum h. We have r\ — (2 — fc)" 1 /! 2 / A. The dimensionless 
parameters L and m are functions of e. Orbits of small ec- 
centricity have m = yk, while those with e = 1 have m = fc. 
We find the orbits in the — V 2 In r potential from the limiting 
case fc — *• 2. 

In Appendix A we show how to improve the accuracy 
of our orbits via perturbation theory; however for most pur- 
poses the simplicity of the initial approximation ^ out- 
weighs the extra complication that accompanies greater ac- 
curacy. Our methods can be applied to non-power law po- 
tentials fsee lKalnaisI (Il979l )) but here, for simplicity, we limit 
ourselves to power laws. 

While our methods can be extended to unbound orbits, 
the results are less pleasing so they are consigned to Ap- 
pendix B. 

In section [3] we use the transformation theory of New- 
ton, iBohlinl [Arnold! and others to transform our orbits for 
< fc < 2 into orbits in potentials with positive powers of 
r. We show how that theory can be extended naturally to 
give a set of transformations that form a group. We develop 
the subgroup of switch transformations and show that or- 
bits in the potentials ip oc r k ~ 2 are conjugate to orbits in 
the potentials r W-k)/k ^-k yk/(2- k ) ^-A,k and r -«/(a-*)_ 

These transformations are not restricted to power laws, al- 



though special simplifications occur for them. Applications 
are made to Plummer's law. 

It is shown that the full group has a transformation that 
connects the Keplerian potential to the isochrone. 



2 ANALYTIC ORBITS 

2.1 General Orbits in Potentials with < fc < 2 

Those looking for orbits in potentials with powers fc outside 
the above range should consult section [3] 

A general orbit of specific energy e and specific angular 
momentum h in the power-law potential ip — Ar k ~ 2 has, in 



the usual notation, r 4> 
r = \/2e + 2Ar k ~ 2 



h and 



h 2 i 



Now 



d<p = (j>dt — hr dr/r ■■ 



dr 



rV2eh- 2 r 2 + 2Ah- 2 r k - 1 



(4) 



In place of r we shall use a dimensionless variable which 
generalises the u (= i), so useful in the Keplerian case: 



u = h 2 /(Ar k ) . 

We also define a dimensionless energy 



E = 



(5) z 



2\ (2-k)/k 



(•») 



(6) 



Now fc dr/r = —du/u so we may rewrite (|4]) in terms of it: 
fc d(j> = - [S(m)]"5 du , (7) 



where 

S(u) = 2Eu + 2u - 



(8) 



and a = 2(fc — l)/fc, which is less than one for < fc < 
2. For the marginally bound orbits e — 0, the u a term in 
((8| disappears so we may integrate (O exactly. Substituting 
u — 1 + cos rj reduces to fc dcj> — dr/. Choosing the zero 
of <j> at that pericentre where r\ = 0, we have r\ = k(j>, so 
the solution for the orbit is of the form of equation Q with 
e — 1 and m = fc: 

u = (£/r) k = I + cos (k<f>) ■ 

Indeed, it was this result that motivated our choice of u as 
the basic variable. Nearly circular orbits can also be nicely 
treated in terms of u, so this encouraged us to conjecture 
that all bound orbits can be found analytically to good ac- 
curacy. 

Our analytic strategy for integrating equation ((7} more 
generally is to replace S(u) with a quadratic function Sq(u), 
which has precisely the same zeros, u a and u p , corresponding 
to the apocentre and pericentre of the orbit. Thus Sq = 
q 2 (u p — u)(u— u a ) whose q 2 , the coefficient of — u 2 in Sq, is 
to be determined so that in some average sense Sq(u) is a 
good approximation to S(u) in the radial range u p ^ u ^ u a 
occupied by the orbit. For example, we find that choosing 
q 2 so that f u " p S(u)u~ 3/2 du = S Q (u)u~ 3/2 du gives a q 
that is good to 2% accuracy. A better choice, given later, 
is the natural starting point for the perturbation theory of 
Appendix A. 

Once S(u) in Q has been replaced by the quadratic 
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Sq(u), the integration is easy. From ((2j), e = (u p — U a )/{u p + 
u a ), so one sets u p = u(l+e) and it follows that u a = 5(1— e) 

O r Q 9 / 2 \ "1 

so Sq = 9 [e m — (« — it )J . If we make the substitution 
u = u(l + e cos 77), we find that the integration of (0 gives 
qk(j> = 77, so the orbit is 

{t/r) k = u/m = 1 + ecos(qk(p) , 

which is of the form {1} with m — qk. 

Crucial to this method of solving for the orbits is the 
knowledge of u v and u a . A critical step in finding them is to 
regard the r v and r a of an orbit as given, in place of its energy 
and angular momentum. Those can easily be found if r p and 
r a are given, but solving the other way around is usually 
difficult. Once the orbit has been found, this approximation 
allows us to determine the radial action and hence the time 
from pericentre to a given point on the orbit. 

Having outlined our general procedure, we now turn to 
solving the circular and nearly circular orbits using u, rather 
than r, as the variable. 



2.2 Nearly Circular Orbits 

For the circular orbits of angular momentum h, the cen- 
trifugal force balances gravity, so h 2 r^ s = —dip/dr = 
(2 — k)Ar k ~ z and so for them u = u c = 2 — k. Also since r 
is zero for them their energy is e c where 

£ c = 5^ 2 / r c ~ Ar k ~ 2 . 

Also S{u c ) = 0, so in our dimensionless variables, c.f. equa- 
tion (El), 



£ c = -ifc(2-fc) (2 " fc)/fc 



(9) 



We consider first orbits with energies not much above 
E c and we set 



(E - E c ) /E c , 



(10) 



then A is small for nearly circular orbits and one at the 
energy of escape. 

E = E c (1 - A) , 

where E c is given by ((9}. For the nearly circular orbits, we 
expand the obstreperous u" term in {§} about u — u c , omit- 
ting terms higher than quadratic in u — u c : 



u — u n 



1 + cr 



Uc 

u — u 
u c 



1 / u - u c 



so inserting this result into |8]l. 

S{u) ~ ku c A + 2A (k — 1) (u — u c ) — q 1 (u — u c ) 2 , 

where again the coefficient of —u 2 in Sq is q 2 and here 

q 2 = [l + A(k~l)]/k. (11) 

Completing the square on (u — u), we have 

S ~ q 2 [e 2 u 2 - (u - u) 2 ] , (12) 

where u — u c + (k — l)A/q 2 and 



e = Aq u [u c kq + A(k — l) z 



(13) 



Integrating with S given by (|12|) by writing u 
u (1 + e cos rf) yields 



kq4> = 77 ; 



(14) 



so the orbits take the form ([T]) with m = kq. Notice that as 
A — ► 0, q -> k~ 1/2 and m -> Vk. 

Whereas these formulae have been derived by neglecting 
the (u — Uc) 3 and higher terms in the expansion of it CT , it 
should be realised that the coefficient of the u a term itself 
vanishes at the energy of escape. Thus, despite this neglect, 
our formulae are exact, not just for A small, but also at 
A = 1. Indeed at A = 1 we see that q 2 = l,u = 1 and 
e = 1. However, even such partial reassurance should not 
deceive us into believing that formulae (|ll[l and (|13|) are 
good enough at intermediate values of e. 



2.3 Analysis of General Orbits 

In the non-linear regime, we see from the e = 1 orbits that 
u has its mean at 1 rather than at 2 — k. It makes little 
sense to expand u a about u c — 2 — k when A is not small. 
Nevertheless we would like to quadrate S, that is, approx- 
imate S(u) by some quadratic function. We adopt a very 
different procedure in the non-linear case. In place of fix- 
ing the energy and the angular momentum of an orbit and 
then determining its shape and size, we choose, instead, a 
pericentric distance r p and an apocentric distance r a - From 
these it is simple to find exactly what energy and angular 
momentum are needed. Equivalently we can fix the values 
of I and e so then r p = £/{l + e) 1/k and r a = £/(l - e) 1/k , 
as can be seen from |T} with mcf> equal to first and then 7T. 
Since r — at both r p and r a , we have for e < 



£ + Ar„ 



- ±h 2 r p 2 = 
e + Ar k a ~ 2 - i/i 2 r~ 2 = 



(15) 
(16) 

which we may solve for h 2 and e in terms of r a and r p or, 
alternatively, in terms of t and e: 



2A 



= i 



k {r P /l) k ~ 2 - {ra/t) k ~ 2 
(r P /£)- 2 - (r a /i)- 2 
fc (l + E )P-*)/*_(l_ E )(»-*)/* 



(l + e y"/*_(l_ E )a/* 



(17) 



Multiplying (|15|l by r p and subtracting from it r„ x (|16(l . we 
deduce 



-f r k — r k 

& 'a ' v 



A 



,h- a (l-e^-q + e)- 1 
(l-e)- 2/fc -(l + e )- 2/fe 



2e(l -e 



2\(2-fc)/fe 



(l + e) 2/k -(l-e) 



2/k 



(18) 



Another alternative, which is the most useful one in the 
equivalent problem in quantum mechanics, is to consider h 
and e as given. Then, eliminating £ in (|18p in favour of h as 
found from (|17[) . we obtain 



A \ A 



2 N (2-fe)/fe 



5(e) ; E = -g(e) , 
where, setting 7 = (2 — k)/k, 



4 D. Lynden-Bell, S. Jin 




Figure 1. The function g(e) for various values of a = 2 — k. g(e) 
is minus the dimensionless energy E. As a — > 0, the graph tends 
to the line g(e) = 1. 



5(e) 



,2/fc 



e(l-e') 



2\T 



[(1 + e) 



[(l + e) 2 / fe 



(19) 



Despite its strange appearance, gr(e) is not a complicated 
function. For k = 1 it is |(1 — e 2 ) and for fc = 2 it is 1. We 
plot (? against 1 — e 2 for several k values in figure [T] 

We wish to approximate S by a quadratic in u which 



must vanish at u 
form 



and u 



so it has to take the 



Sq 



q (u — u a ) (u p — u) 



2 r 2-2 

q e u 



(u- 



(20) 



where q is yet to be determined and has been given that 
notation to conform with our earlier definition that — q 2 is 
the coefficient of u 2 in 5*. In the above, 



I , , h 2 

U = - (U a + Up) = 



2A \ p 



k +r a k 



(21) 



from which we see that u is twice the final expres- 
sion in (|17|l but without the £ . The 'eccentricity' e is 
(u p — u a ) I (tip + u a ) as always. Using (120[) for S(u) with the 
substitution u = u(l + ecos77), we readily integrate equa- 
tion to obtain qk(j> — r\ as before. So the solution is still 
equation {1} with m — qk, but we must still determine q 2 . 

A useful approximate formula, good to about 2%, is 
given by setting S(u)u~ 3/2 du = S Q (u)u~ 3/2 du. 
This gives, for a 1/2 (i.e. a / 2/3), 

, i(2- 



(2-(a-i)- 1 ) + ((a-i)- 



)(1 + iyr^2) 



1- V^e 2 ) 



(22) 

with m = qk as before and u can be expressed as a func- 
tion of k and e only, via (|17|l and (|21[) . We have chosen this 
power of u in the integrals we equate above, as it gives the 
best agreement to the true m without compromising on the 
simplicity of the expression for q. We note that choosing it -1 
in the integrals gives as good an agreement as using u~ 3 ^ 2 , 
but with the resulting m becoming overestimates on the true 
value for a < 1 and vice versa for a > 1. Choosing an expo- 
nent between —1 and —1.5 results in better agreement still, 
but we lose the simplicity of the resulting analytic expres- 
sion for q. The exponent of —4/3 is as good as any, but gives 
the following somewhat awkward result: 

2 5(<t — l)(e + + e 

q = 



+ (2 -a) [2(e + + e_) 



(3(7- l)(e+ + e- 



where e± = (1 ± e) 1/3 . 

The angle between successive apocentres is important 
as such angles accumulate as the orbit is prolonged. We now 
determine q to get this angle as accurately as possible. It is 
given by 



* 2tt 2 

$ = — = T 

m k 



s 



-1/2 



du = 



2 p (Sq 

k Ju a V s 



1/2 



du 



S, 



1/2 



1 

m 



Sq 
S 



1/2 



dr) , 



hence the average of (Sq/S) 1 ^ 2 over r/ must be one. We 
evaluate this average over eight points around — n < r\ < n. 
There is a difficulty in evaluating (Sq/S) 1 ' 2 exactly at the 
apocentre when e = 1 since the apocentre is at infinity. 
Surprisingly, the result of taking the limit of r a as e — > 1 
gives a different (and wrong) result from setting e = 1 and 
then evaluating the limit as cos rj — > — 1. We get around this 
by using cos?; = —0.990, where everything is finite, in place 
of r\ — ty. 

At pericentre u p , both Sq and S are zero but the limit 
of USq/S) 1 ' 2 is 



2ue 



cr(2 — Up) — 2(u p — 1) 
We evaluate 
1 



(u p — u)(u ■ 



2Eu" + 2u 



where u = u(l + e cos rj), at the other seven points ±tt/4, 
±7t/2, ±37t/4 and cos 77 = —0.990. For convenience, we label 
the values of u at ±7r/4 and ±37r/4 as u = w(l+e/v2) = u+ 
and u — u(l— e/V2) = u~ respectively. Our estimate of 1/q 
is the average over the eight values that result: 



8 ^ 



8 r, w Ml/2 

(Up — Ui)(Ui — u a ) 



S(ui) 



kq 



(23) 



At each a, the resulting q is a (somewhat complicated) func- 
tion of e, since E, u a , u p , u axe all functions of e. 



2.4 Comparisons with Computed Orbits 

Orbits were computed in the x,y plane from the Cartesian 
form of the equations of motion for potentials with a = 
0.25, 0.55, 0.75, 1.5 and 1.0. The last provides a valuable 
check that we get m = 1 in the Newtonian case, even at 
very high eccentricities of order 0.999. We also checked that 
m — y/k for nearly circular orbits and that m — > k as e — > 
1 for all the values of a. As m varies quite rapidly with 
eccentricity as e 1, accurate computations are required 
at high eccentricities. Orbits in the logarithmic potential 
which gives a constant circular velocity are considered later, 
in section [231 the approximation adopted there is somewhat 
different. 

Figure [5] shows a comparison between our estimated 
values of m and the computed values of m for potentials 
with a = 0.25, 0.55, 0.75 and 1.5. For all values of a shown, 
the deviation of m calculated from the analytical formula 
(|23[) is only a fraction of a percent. Notice that both plots 
are against yl — e 2 so that high eccentricities are on the 
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0.5 



-0.5 



1.5 



0.5 



-❖a 



a a = 0.25 

□ a = 0.55 

* a-0.75 

x a=1.5 



- -A - . 



<x-X' 




a=0.75 



m true 



0.2 



0.4 0.6 

(|- c 2)l/2 



0.8 



Figure 2. The top panel shows the percentage difference between 
the true value of m from computed orbits and 77123 calculated 
from the analytical formula {23} . plotted as functions of vl — e 2 
for different values of 2 — fc = a. The largest deviations occur at 
high eccentricities; the region 0.1 to 0.3 in Vl — e 2 corresponds 
to 0.995 > e > 0.954. The bottom panel shows the values of m, 
estimated using (1221 and J23}, along with the true values, again 
for four values of a. 



left and low eccentricities are on the right. It is, of course, 
possible to read off m as a function of y/l — e 2 or of e from 
the computed points in this figure. 

Panel a) of figure [3] shows a computed orbit in the po- 
tential with a = 0.25 together with an orbit of the same 
m, r a and r p but calculated from the equation {t/r) k = 
1 + ecos(mcj>). This demonstrates how the shape given by 
equation Q fits the computed orbit. A better fit is obtained 
using the perturbation theory of Appendix A. The dotted 
orbit in panel b) is (£/r) k = 1 + e cos(rn[^ffii>) and the gradual 
precession due to the error in the estimated 77 ^23] is readily 
seen. 

Figure |3] shows two orbits with the same ratio of r a /r v 
but in the potentials with a — 0.55 and a — 0.75. Because 
the definition of 'eccentricity' we gave in equation ([2]) de- 
pends on a (through fc), these orbits have eccentricities of 
0.662 and 0.596 respectively. Notice that the two drawings 
have the same number of apsides but these have precessed 



2 - 



>> 



-2 




Figure 3. An orbit in the 2 — k = a = 0.25 potential, chosen to 
be in the high eccentricity region where the approximation is least 
good. The full line is the computed true orbit. In panel a), it is 
compared to (l/r) k = 1 + ecos(mif>) (dotted line) with m chosen 
to agree with the computed orbit. This shows the deviation in the 
shape of the lobes. In panel b), the value of m is estimated by 
the 8-point average of equation l|23j l and the error in m is readily 
seen from the different precession of the solid and dotted orbits. 
The difference in the true value of m and that estimated from the 
8-point average is less than 0.5%. All orbits start at the starred 
location. 



-2 







/ \ a = 0.55 


a — 0.75 











-2 



Figure 4. Two computed orbits with the Same ^min/^max in 
potentials with different 2 — k = a are compared with the analytic 
orbits (£/r) k = 1 + e cos(m tmo (f>) (dotted lines). The a = 0.55 
orbit has e = 0.662 while the a = 0.75 orbit has e = 0.596 and 
precesses much less rapidly. Both the shapes and the precession 
rates of these orbits are well represented by the analytic formula. 



much less for the a = 0.75 orbit as the potential is closer to 
the Keplerian a = 1. 

Figure [5] shows two orbits in the a — 1.5 potential for 
which the precession is forwards because a > 1 whereas the 
other illustrations all have a backward precession. Under the 
transformation £ = z k ^ 2 considered in section|3] these orbits 
transform into ones in the potentials tp oc r 2a ^ k = r 6 . For 
orbits with a > 1, the straight perturbation theory giving 
m = kq with q given by ((23} yields m to better than 0.5%. 

It is often useful to have a vectorial way of delineat- 
ing orbits and the velocities of particles describing them. 
To do this, we generalise Hamilton's eccentricity vector 
which has magnitude e and points toward pericentre. As 
our pericentres precess within the orbital plane, we in- 
vent a rotating eccentricity vector. If we start at pericen- 
tre with e = eo we take e at later times to be given by 
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>> 




Figure 5. Two orbits in the 2 — k = a = 1.5 potential. The orbit 
on the left has a forward precession by nearly 180° whilst that 
on the right has a forward precession by a little over 270°. The 
turn close to the origin in the latter cannot be seen but is a more 
rapidly turning version of that seen on the left. The lobes are 
numbered in sequential order for the higher eccentricity orbit. 



e = eo cos [(1 — m) <f)] + h x eo sin [(1 — m) <j>\. This e obeys 
de/d(j> = (1 — m) h x e. The angle between the radius vector 
to the particle and the eccentricity vector is then m<j> and 
the equation of the orbit |l} can be rewritten 



(£/r) 



1 



(24) 



The transverse velocity of the particle is clearly hxf/r and 
the radial velocity can be obtained from the orbit and the 
energy equation. Using our approximations quadrating the 
latter, we find 



4 [hxf+ q(r/£) k h. (e x f ) f J 



(25) 



Given v and r at one time one might wish to use these 
equations at a later time. Then one needs to find e, h, £, q 
and m from the initial v and r together with the known 
potential ip = Ar~ a . From v and r it is easy to con- 
struct e = t;« 2 — Ar~ a and h = r x v, from these E is 
found. For given E and a, e may be found from figure Q] 
£ then follows from (117[l and m,q from (|23[) . The direction 
of e within the plane perpendicular to h then follows from 
e.f = {£/r) k — 1 with the ambiguity in angle resolved from 
v.f = r~ 1 q (r/£) k h. (e x f ) which follows from the v equa- 
tion above. Thus all the orbital parameters are determined. 



2.5 Action, Adiabatic Invariants and Time 

So far we have concentrated on the shape of the orbit in 
space, however the time from pericentre to any point of the 
orbit is just as important. Both can be obtained from the 
action function S r , whose relationship to S(u) is given be- 
low. 



S r = I p r dr — J fdr 

= J \/2e + 2Ar-°> - h?r~ 2 dr 



h Lf2eh~ 2 r 2 + 2Ah- 2 r k - 1 r' x dr 



k^h / ^n" du . 



(26) 



If we now use our quadratic approximation we find 
S r = qk~ 1 h / \J e 2 u' 2 — (u — u) 2 u~ 2 du ; 

J u 

setting u = u (1 + e cos 77), this becomes, setting / = 1/e, 

1 f v sin 2 7/ 

S r = qk h — 3 dr] . 

Jo (/ + cos ry) 

Now the related integral 

(1 - cos 2 n ) d v = fv i-f 

f + cosr) J V/ + cos?i 



+ / — cos 77 ) dr] 



-2^/p - ltan" 1 



f + 1 \2 



+fV - sin 77 , 

and the integral that we want is just —d/df of this, so 



/ 

■j 



v sin 2 77 



(/ + cos 77) 



-dr] 



2/ 



: tan 



tan - 



/ + 1 



sm 77 



/ + cos 77 

so, putting this in S r and remembering that / = 1/e, 



2 
(27) 



S r — qk h 



-r) + 



tan 



1 — e 77 
tan — 

1 + e 2 



+ - 



e SU177 



1 + e cos 77 



The adiabatic invariant is given by 



(28) 



(29) 



Now 

dJ r /de\ h = ±-<£±dr = P r /{2ir) , 

Z7T J r 

where P r is the radial period while 

-dJr/dhL = -!- ihr~ 2 \dr = ^- * <pdt = */(2tt) . 
2tt J r 2tt J 

J r /h is a function of eccentricity, so the partial differentia- 
tion is best done via 

d{J r /h)/de = d(J r /h)/de (de/de) h 

and 

(8Jr/dh) e = Jr/h + A (J r /h) h f^- 



(dS r /de) h = J± 



= t 



(30) 



so this expression gives the time to any chosen point in the 
orbit. In practice, S r /h is a function of e and 77 so the partial 
derivative is done using 

(dS r /de) h = (dS r /de) h {de/de) h . 

In the general case, use of — oc r k as a variable does 
not lead to a prettier equation for t, such as the one Kepler 
derived for k = 1, but see the next section for logarithmic 
potentials. 
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In the equatorial plane the total action is Sa = S r + 
h z (j). The action variables are J r and h z — h. The angle 
variables are the phases of the oscillations in r and <j> and 
are given by w r = dSA/9J r and = dSA/dh. For general 
orbits, the action variables are most often employed when 
the potential is of the more general separable form; ip = 
Ar~ a — B (6) /r 2 ; h 2 is no longer conserved but I = h 2 — 
2B (6) is. The general action is then Sa = S r + Se + h z <f> with 
Sg = /* yi - h 2 cosec 2 6 d(9 and J 9 = ^ § dSg/dd dQ. The 
angle variables are we = dSA/dJg and io^> = dSA/dh z . 



2.6 Logarithmic Potentials a — ► 

For small a we write ip = Ar" a = Ar^e"" ln(r/ro) and 
expand to obtain ip = Arp [l — a In (r/r ) + (a 2 )] • We 
set >1 = V 2 /a and consider taking the limit as a — > while 
keeping V 2 fixed so A tends to infinity. To keep a finite 
potential, we have to subtract the constant Ar^ 01 from ip, 
so we obtain a new potential 

* = ip - Aro a = -V 2 In (r/r ) . 

To apply the methods of section |2~T1 we define 

u = h 2 / (V 2 r 2 ) 

and consider orbits defined by pericentric and apocentric 
distances r p and r a - In place of equations (|15|l and (|16l) we 
then have 



V 



M(r a /r v ) _ 1 2 2 



1 + e 



(31) 



so itp = (1 + e)w , where u — i In [ - — 6 ) , and (32) 

2e \ 1 — e . 



V' 2 



2 r 2 In (r Q /r ) - r 2 In (r p /r ) 



[(1 - e) In (1 + e ) - (1 + e) In (1 - e) + 2eln (^) 



(33) 



The orbital equation reads 
dr 



d<f> = 



~du 



ry/2eh- 2 r 2 - 2V 2 h~ 2 r 2 In (r/r ) - 1 2^Sdu) 



where 



Sl {u) = 2Eu — u In (mo /u) — u 2 , 

= ^[(l-e)ln(l + e )-(l + e)ln(l-e)] 



+ «ln(|)- M 2 (34) 



and E = e/l/ 2 which is given in terms of e via (|33[1 . We now 
approximate Sl(u) by the quadratic (|20|) . which shares the 
same zeros. 

As in section 12.31 a useful analytic expression for q is 
given by again equating 



/ S L (u)i 



-' i/2 du=Au 1/2 (^/TT7 - VT~ e) 



[|(2+ X /T _ 



1 .9 FF 1 " 



1.8 T 

« m. 



1.6 - 



1.4 



0.2 OA 0.6 0.1 

(, _ e 2)i/a 



□3D 
1 




0.2 0.4 0.6 0.8 1 



Figure 6. The left panel shows a comparison between the true 
value of m (open circles) and that estimated from equation (135 D . 
which gives m to better than 2%, and the 8-point average in 
equation (1406 . which gives m to better than 1%. The right panel 
shows (in thick solid) the true value of 1/q = k/m overlaid (in 
crosses) with that derived from the 8-point average. The agree- 
ment is clearly good to a fraction of a percent. The individual 
contributions made by each of the points given in the terms of 
equations l|39|l and l|40j l are also shown. From the top, the vari- 
ous lines correspond to u a (dotted), (dot short-dash), u (long 
dash), u- (short dash) and u p (dot long-dash). 



to 



S Q {u)u- 3/2 du=^q 2 u 3/2 (vTTi- VT~^ e) x 



(l- \/l-e 2 ) 



where u is given by (|32p . This yields 
2 + Vl - e 2 - 3/u 



2 

q = 



2 (l - vT^i 



ra = 2q 



(35) 



This expression is good to 2%. Once again, more accurately 
we calculate g _1 from the average of Sq / Sl over 77, where 
u = u(l + ecosrj). At u p , 

2E = ln(uo/u p ) + u p ; 
dSL/du\p = 2E — ln(lto/u p ) + 1 — 2u p = 1 — u v ; 

dSQ/du\ p — ~2q 2 ue ; 



so \/ Sq/Sl — > q\/2ue/(up — 1) = g/^ p 

Ip 

and at 77 ± 7r/2 , \/Sq/Sl = queS^ 1 ^ 2 , 



(36) 
(37) 



where S 1 = Sl[u). As found in section T2.3I the apocentre 
once again poses a problem and we evaluate 



Ha = 



(Up — u)(u — Ua) 



Sl(u) 



(38) 



at u = u(l + Xe) where A = —0.990, as before, and where 
Sl{u) is given by (|34|l . So a 4-point estimate of q is given 
by 



q 1 ~ g 4 1 = \fl p + fia + 2ueS 1/2 j /4 , 
and a 8-point estimate is likewise 



«8 



x = [4g 4 - 1 + v / 2^(s; 1/2 + 5: 1/2 ; 



(39) 



(40) 
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Figure 7. An orbit in the logarithmic potential. The solid line is 
the computed orbit, and the dotted line shows the approximate 
orbit (i/r) 2 = 1 + ecos(m tme 0). Overlaid (and hardly visible) is 
a dashed line showing an orbit with the estimated m, calculated 
from l|40jl . instead of mt rue - The error in this value of m is less 
than 0.1%. The heavy line shows the piece of such an orbit for the 
trailing Magellanic Stream, with the location of the Magellanic 
Clouds indicated by the open star. 



where S± = Sl(u±) and u± — u(l±e/v2). Figure [6] shows 
the contribution of each point in (|39[) and (|40p along with 
a comparison of the true m and the value derived from the 
8-point estimate via 7T f4Q] = qsk. Figure [7] compares an ap- 
proximate orbit to a computed one. 

The time to a given point in the orbit is given by 



t= f r *: : 



h 



'2qV 2 
h 

' 2uqV 2 
h 



-du 



uy (ue 



(«• 



e cos t) 



tan 



uqV 2 yr^2 

so the radial period is given by 



1 + e 



tan (q<f>) 



(41) 



(42) 



However, a much more interesting result comes from fol- 
lowing Kepler, whose equation comes, not from integrating 
the u equation directly, but by first making the substitution 
v = 1/u = V 2 h~ 2 r 2 . This gives 



dv 



2qV 2 uVl - e 2 J Vp ^/v 2 e 2 - (v - v 2 ) 2qV 2 uVT^e 2 




i 



time [arb. units] 



Figure 8. The almost simple harmonic variation of r 2 as a func- 
tion of time in an orbit with r a /r p = 2 in the logarithmic poten- 
tial. The dotted line is a sinusoid of the same period. 



where we have written v 
setting k = 2qV 2 u\/l — 



= v(l—e cos x) and v 
e 2 hT x we see that 



X — Kt and r 



h 2 V~ 



1(1 



[1 



3(/S*)] 



= S(l-e 2 ); 



(43) 



so with this approximation r vibrates harmonically. Fig- 
ure[S]shows the computed r 2 (t) for an orbit together with the 
harmonic approximation. The adiabatic invariant is given by 
J r — § rdr ~ § \J S(u)u~ 2 du. This integral was eval- 
uated in equation (|28p . so for this case k = 2 and 



-hq 



1 



(44) 



It should be emphasised that while we have set ourselves 
the target of getting analytical formulae that give m to 1% 
or better for all eccentricities, we have not paid attention to 
minimising errors in the temporal periods. We find that such 
errors are indeed higher and no doubt our formulae could be 
improved upon by a study of such errors. 



3 TRANSFORMATION THEORY 

Newton (1687) realised that the ellipse was a possible orbit 
both in a harmonic central potential and in an inverse square 
law. In the first case the centre of force is at the centre 
of the ellipse, while in the latter case it is at the focus. 
This led him to pose the question under what circumstances 
can the same curve be the trajectory of a particle unde r 
a force from one of two different centres. Newton's 
discussion of this is well d escribed in Ch a ndrasekhar's 
book, as is later work bv lBohlml ( 191ll ). lLevi-Ci vital 
and l Arnold! (|l990l ). All demonstrate the transformation that 



1714) 



1995) 



19241 ) 



conve rts the harmo nic ellipse into the Kepler ellipse and vice 
versa. ICollasI (|l98ll ) gave the relationship between equivalent 
potentials. Here we show that this transformation, Si, can 
be embedded into a larger set of transformations that form 
a group. We mainly concentrate on the subgroup of switch 
transformations which have six members but which give just 



three related potentials 



, r and r in the Kepler case, 



since r is self conjugate under one of the transformations. 
In the comple te group, the se potentials are also related to 



in the complete group, the; 
the isochrone (|Henonlll959h 
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The energy equation of a central orbit of angular mo- 
mentum h in a potential ip(r) can be written 



1.2 i . 1,2-2 

-r = ip + e- -h r 



(45) 



Now consider the transformation r = p(r) , dt — 
dt/r(r). Setting F — p V we find 



For this to be an orbital equation like (|45[) , the 3 terms 
on the right must be , £ and — |nf -2 , but not neces- 
sarily in that order. 

The transformation Si that leads to the Newton-Levi- 
Civita- Arnold result is found by taking h 2 — h 2 , identifying 
the last terms but switching the roles of the other two. Thus 
for Si we set 



F 2 ip(r) = i , F 2 e = 4> , F 2 r~ 2 = r~ 2 
so we obtain the transformation 



(47) 



r = r*i>(r)/i , (48) 

with i/>(f) = er 2 /f 2 = ee/ip. Applying this to the power- 
law potential tp — Ar~ a yields r = r 1- "' 2 ^ A/s so ^(r) = 

? 2a/(2- a ) ^ (gr / / A )4/(2- Q )j ^ ote . a = 1 giveg ^ K ~2) ^ 

quantity in square brackets is, of course, constant. However 
we could alternatively leave the first term on the right of 
(|46[) . identifying it with ip but switching the roles of the 
other two terms. This leads us to the transformation S2 in 
which 

1 

2 



F 2 ip 



F 2 e. 



1 1,2 — 2 

— n r 
2 



-F 2 h 2 r~ 2 



from which we deduce ipr oc ip and f 
2 1 hh 



where 



50 S2 is an inversion accompanied by the change in potential 
ip oc f _2 i/> (a 2 /f) oc f a ~ 2 (for power laws). More generally 

we may ask that i (^f) 2 = F 2 \<' 2 but that the three terms 
ip, i and — \h 2 f ~ 2 that constitute this quantity are indepen- 
dent linear combinations of F 2 ip, F 2 s and F 2 (—^h 2 r~ 2 ). If 
one applies two of this more general class of transformations 
one after the other, it is simple to see that the net result is 
a transformation of this class, that the identity transforma- 
tion belongs to the class and that every transformation has 
a unique inverse in the class. These transformations form a 
group in the sense of group theory, since they clearly obey 
the associative law 

(T3T2) Ti = T 3 (T 2 Ti) = T3T2T1 . 

Each transformation now gives a new f (i) in the new 
potential ip(r) that corresponds to the old r(t) in the old po- 
tential ij). To get the transformation of (f> we remember that 
p' — df/dr and for Si, h — h so under the transformation 

51 using (gSJ) and (gTJ) 



dm = — — 

r 2 dr/dt 

but from ESI) i 



h d In f d In r 



din r 



dcj) , 



din r / 1 din ib 

= H - 

din r \ 2 din r 



This equation takes a particularly simple form when f is a 
power of r, for then <f) oc cj>. Furthermore, this occurs if and 
only if ij) follows a power law in r; so 



ip=Ar 



= 1 



V>=£ ( -J ) r a 



1 I 



1 /> 



Notice that the new potential depends on the energy of 
the old orbit, so a pair of orbits of different energies in 
the old potential will map into a pair of orbits in two 
diffe rent new potential s that differ by a constant factor, 
c.f. iRoscmist fc Pucaccol (|l995h . If we write z = and 
z = f e l( * then, from the above, the Si mapping is of the form 
z oc z 1_a,/2 , i.e. a conformal map in the complex plane. In 
general, a closed orbit will map into an unclosed Lissajoux 
rosette, but when 1 — a/2 = N1/N2 where Ni and N2 are 
relatively prime integers, then the transformation of an orbit 
that closes after one turn will be an orbit that closes after 
N2 turns which has Ni times as many apsides. 

For q = 1, we have the famous example that transforms 
Kepler's ellipse into the simple harmonic oscillator. This is 



-2 we find 



r = r = (A/e)* , = \tj> , i> = e (i/A) 2 f 2 
If we apply Si again, this time starting with a 

Si [f] ocf 2 oc r 
and Si ocr -1 

so apart from a possible rescaling, the double transformation 
S 2 leads us back to the beginning. This is true generally, not 
just for power laws, since from equation (|48|) . 



51 [f] oc r\Nj oc r^/v^ oc r and Si oc -L oc tp . 

We shall ignore the dull rescalings in what follows and write 

5 2 = I, the identity. This is in agreement with the concept 
that a repeated switch leads to no transformation. 

We now apply the transformation Si to orbits in one of 
our potentials ip = Ar~ a with < a < 2. The new potential 
will be ■tj) oc f 2a / k , which will be a positive power of f and 
the transformed orbit takes the form 



1 + e cos 



2m - 
IT * 



which is indeed an ellipse when m = fc as for the Kepler case, 
which transforms to the harmonic potential. Remarkably, it 
is always f ~ 2 on the left whatever a we start from, but the 
values of 2m/fc vary with a. 

Si is the basis for the regularization of the close en- 
counters of two bodies c arried out in three dimensions by 
Kustaanheimo & Stiefel (|l965h . 



3.1 The Switch Subgroup 

If we try to find a transformation that switches the angular 
momentum and potential terms in ((46} while leaving the 
energy term unchanged, we fail because e = const = e and 
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with F constant we are unable to accomplish the desired 
switch. 

However, we may apply first Si and then S2: 



Si [tP] = ee/i> ; 



S 2 [Si[r]] = - 



hS 2 



['•] 



e 1 



h 



S 2 [h] 



S 2 [Si 



2S 2 [e]e 2 

. r 

h 2 



This transformation is not the one we obtain by applying 
5*2 first and then Si: 

r = S 2 [r] = \l% 1 -; S 2 [tp] = -ff^r 2 ; 
2 Wee r h 2 



Si [S 2 [r}} = - 



T \fi> 



Si [32 



{-2eSi [e])* 



= -i/i 2 Si [e]/(^r 2 ) 



Applying this double transformation twice gives 
r4 oc 1/ ^r\/^ and tp oc r 2 , 

which is the same as S 2 Si up to constants of proportionality, 
while a further application of Si £2 gives 



re oc r 



and 



ipe oc tp 



so the triple application of S1S2 gives a multiple of the iden- 
tity. 

Before going any further, let's see where we can get if 
we start with ip oc 1/r. We can get to ip oc r 2 and back 
using Si, but using S 2 leaves ip oc 1/r, so Newton's law 
is invariant under S 2 . However, S 2 acting on r 2 leads to 
tp oc f~ 4 , but a further application of Si leaves r~ 4 in- 
variant. Thus under the transformations considered so far, 
there are conjugate orbits in the r _1 , r 2 and r -4 poten- 
tials. More generally, if we start with tp oc r~ a , then Si gets 
us to ip oc r 2a ^ 2 ~ a ' > (where we have dropped the tildes), 
while S2 brings us to tp oc r a ~ 2 . The double transforma- 
tions S2S1 and S1S2 give tp oc r - 4/(2 - Q) and r 2(2 - Q)/Q 
respectively, while S1S2S1 leads to tp oc r~ 4 ^ a . Further 
applications only bring us back to potentials already in- 
cluded. In fact, there are six transformations in this sub- 
group, yielding a conjugacy of orbits in the six potentials 

r -a r 2a/(2-a) ^-2 ,.-4/(2-0,) ,.2(2-a)/a and p -4/a p or 



q = —1, these are r, r 2 ^ 3 , 



-4/3 



and r . For 



a = 1/2, they are r~ 1/2 , r 2/3 , r _3/2 , r~ 8/3 , r 6 and r -8 . 
These powers become somewhat bizarre for small a. 

1 la ■ -1/6 2/11 -11/6 -24/11 22 , -24 

a= 1/6 gives r 1 , r ' , r , r , r and r ; 

i , e • 1/6 -2/13 -13/6 -24/13 -26 j 24 

a — — l/b gives r , r ' , r , r , r and r ; 

degeneracies similar to those for the Kepler potential occur 
for a = l, 4 or ±2. 

The simple relationship <p — <p{l — a/ '2) holds only for 
power law potentials under the Si transformation. Under S2 
we find 

d(p = hf 2 (dr I dt) df — —(—2e) 1 ^ 2 r^ 1 dr/r 
(~2 £ ) 1/2 

= -- — r — rd( t ) > 



which no longer gives a simple relationship of cp to (p. How- 



fr 



^ v ! ^1/2 J 1 "V — X l r J — h l/2 

if x( r ) is known for the first orbit then, with f(r) known, 
x(r) is known for the second. 

These transformations are not restricted to power- 
law potentials. Under Si, Plummer's potential tp = 
H (r 2 + 6 2 ) 1/>2 transforms into 

2eb 2 f- 2 



1/4 



.1/2, 



= ~X, so 



1 ± 



4^ 2 b 2 



while under S2 it becomes 



•0 OC 



(r 2 + b 2 ) 



1/2 



/4ea 2 



^r 2 + 1 



1/2 



V h 2 h 2 

It would be tedious to give the complete set but there are 

six;</>, Si[tp], S 2 [tp], S 2 [Si[tP]], Si[S 2 [tP]] and Si[S 2 [Si[tP]]]. 

3.2 The Larger Group 

When we ask that F 2 r 2 = (df/dt) 2 but in place of merely 
switching the terms on the right of (|46|l we ask that those 
terms are linear combinations of tpe and — |/i 2 r _2 , we obtain 
the full group of transformations. A general transformation 
of the group is then 

tp = F 2 (antp + ai 2 e - ai3±h 2 r~ 2 ) , 

e = F 2 (a 2 itp + a 22 e - a 23 ±h 2 r~ 2 ) , 

-±h 2 r~ 2 =F 2 (a S itjj + a S2 e - a 33 ±h 2 r~ 2 ) , 

where J2 n a nj = 1 for j = 1, 2, 3. 

If we take the particular transformation with 1231 = 
032 = ai3 = «23 = «12 = then 021 = 1 — an and 
a 22 = 033 = 1 so we get, taking h = h without loss of 
generality, 

tp = F 2 autP , e = F 2 [(1 - a u )V> + e] , r~ 2 = F 2 r~ 2 . 

The second of these gives the relationship of f to r when the 
value of F 2 is taken from the third: 

f 2 = r 2 [(1 - an) tp + e] /e . 

Taking tp = 
to find r(f] 



Taking tp = GM/r as our initial potential, we readily solve 
if 2 (Vr 



+ b 2 -b), 



where b - 



(l-an)GM 



' ee 



GM 



The potential tp = ail< ^^ T = an ( — ] , 

V r 2 \ej Jf 2 + b 2 + b 



where we used f 2 = (V f 2 + b 2 — b) (\/r 2 + b 2 + b) . 

The potential tp is the isochrone, see iHenonl (| 1959h . 
This is the most general potential in which all orbits can be 
found using only elementary functions (trig o nome tric etc.) 
as stated bv lEggen. Lynden-Bell fc Sandagei l|l962h ; the de- 
tailed proof of this was only published many years later in 
lEvans et all (|l99Ch . 



4 CONCLUSIONS 

We have found crude, but useful, approximations to Abelian 
functions by quadrating the expression under the surd while 



Analytic Central Orbits and their Transformation Group 11 



keeping the end points as the constant parameters. For our 
problem, these methods give accuracies to better than 1%. 

We have shown that the e = 1 'parabolic' orbits at the 
energy of escape can be solved exactly, and we have given 
analytic expressions for m(e) which hold for all eccentricities 
^ e ^ 1. We have thus illuminated why Struck found these 
orbits to be such good approximations at low and moderate 
eccentricities. 

The transformation theory has allowed us to extend 
these results to orbits in potentials which are positive pow- 
ers of r and we have extended the transformations to form 
a group. 



«3 = | [p(0) -p(tt) - v/2(p(7r/4) -p(3ir/4))] ; (A12) 
«4 = l[p(0)+p(7r)+2p(7r/4)] . (A13) 
The perturbed orbit is given by the implicit equations 

4 

m<j) = ri + '^^n~ 1 a n sm(nr]) and u = u(l + e cos 77) . (A14) 

n=l 

In the above solution for the a n , we have not used the combi- 
nation i [p(7r/4) +p(37r/4)] because the condition < p >= 
ensures its consistency and we replace p(n) by p(<f>) with 
<j} = cos _1 (-0.990). 
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APPENDIX A: PERTURBATION THEORY 

We have S(u) = 2Eu a + 2u - u 2 ; a = 2(1 - a)/(2 - a), 
< a < 2. We rewrite it in the forms; 



S(u) = q 2 [(ue) 2 -(u-u) 2 ]/[l+p(u)] 2 
= q 2 (u p — u) [u — u a ) I [1 + p(u)] 2 , 



(Al) 



where p(u) is the perturbation function that allows for the 
difference between S(u) and our quadratic approximation to 
it. The orbit is given by 



qk d(j>=q 



-du 



(l+p(u))du 



(u — it) 2 



(1 +p[u(l + ecos??)]) dr) , 



where u — u (1 + e cos 77). 

Now p may be expanded in a Fourier series in 77: 



(A2) 
(A3) 



p — ao + ai cos 77 + 02 cos 277 + 03 cos 377 + 0,4, cos 477 + . . . (A4) 

where a n = A cos (77,77) dr/ for 77 7^ 0. We chose q so that 
the average of \J Sq/S over 77 is 1. This ensures that the 
average < p >— 0, so ao = 0. 

Keeping just those terms with n 4 we find that ^(77) 
is given by 



p(0) = ai + a 2 + a 3 + a 4 

p(7r) = — 01 + (72 — 0.3 + (74 
p(7r/2) =— a2 + (74 

p(7r/4)=oi/-v / 2- 03/V2- a 4 
p(3 7 r/4) = -ai/v / 2 + a s /V2 - a 4 



(A5) 
(A6) 
(A7) 
(A8) 
(A9) 



From these, we may deduce the following: 

ai = 4 [p(0) - p(tt) + V2 ( p (tt/4) - p(3»r/4))] ; (A10) 



a2 = i[p(0)+p(7r)-2p(7r/2)] 



(All) 



APPENDIX B: UNBOUND ORBITS 

When e > the E term dominates at large distances (u 
small). Indeed when E > 1 the potential term never domi- 
nates. We write the equation for the orbit in the alternative 
forms, 



k d(j> - 



du 



du 



^j2Eu a + 2i 



N /£ fc / Q (7 2 ( 1 - Q )E(f/) 
kE k/{2a) dU 



where E (U) = 2 + 2U a - E k/a U 2 



(Bl) 



and 



U = E 



-l/ct 1/fc 



E- 1/a (h'/Af'r- 1 



When E > 1 we use the U form everywhere and ap- 
proximate the U" term in E as a quadratic. We specify our 
orbit by the values of the impact parameter b = h/\/2e and 
the value of the perihelion distance r p . The energy equation 
at r v is 



e = -Ar p a + \h 2 r~ 2 

hence 1 + Ar p 2 e~^ = b 2 /r 2 , = f3 2 



say; 



b and r p are both specified and A is known, so we find e as 
e = Ar; a / (/3 2 - l) , 



and with e now known h is given by h = 6y2e. 

From (|6]) we can now deduce the dimensionless energy 



/ t \ otlk 



[(V2/))"/(/>'-l)] 



We may also express U in dimensionless combinations 
U=(E)~ 1/a [2p 2 /(p 2 -l)] 1/k (r p /r) . 

We require our quadratic approximation to E to be ex- 
act at r p , that is, U = U p and exact at U = and at the 
centre of the range U p /2. Then the approximation takes the 
form 

E=2 + 2U a - E k/a U 2 ~ {2 /U P + Q 2 U) (U p - U) (B2) 
= Q 2 [e 2 U 2 - ((7-f7) 2 ] , (B3) 
where 

Q 2 = (2/U P f +2(2/U P ) k - E k/a ■ U= iU p - (Q 2 C/ P ) _1 ; 
e 2 = l + 2(Q 2 U)~ 1 . 
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Thus for E > 1 we may integrate, using this quadratic 
approximation to obtain " = - = 1 + e*cosm,0, where 

ml = Q 2 E~ k / a and i = E' 1 ' 2 [2/3 2 / (/3 2 - l)] 1/k (r p /U) = 



When < E < 1, the _B term dominates at large r 
(small u) but the 2u term dominates it when 1 < U. We ap- 
proximate the smaller of these terms in each region and make 
sure that the two approximations to S(u) join smoothly with 
the same gradient at U = l,u = E k ^ a . The pericentre lies 
in the U > 1 region where the S(u) form is appropriate, so 
we need S(u p ) = 0. 

S(u) ~ (co + ciu)(up — u) = c\ [e 2 ?t 2 — (u — u) 2 ] ; u i5 fc / c 

where Co and ci are constants to be determined and e, u 
follow them. 

Since U p no longer lies in the zone where the E form is 
used, our former approximate form (|B2[) for E is not appro- 
priate. We write, instead 

E(t/) = 2 + CiE/ - CW 2 = C* 2 [e 2 t7 - (U - I7) 2 ] , 

where Ci and C2 are constants. 

Now looking at JUJl, it is S(u) and £ fe/cl ?7 2(1 - Q,) E ((7) 
that have to be continuous with a continuous derivative at 
the junction point u = E k ^ a , where we demand the deriva- 
tives be exact. We have du/dU = kE k ^ a there. Continuity 
requires: 

(co + ci£ fc/a ) (u p - E k/a ) = E k/a (2 + Ci - C 2 ) , (B4) 

whilst the condition on the derivative provides: 

(a + 1 - £ fe/a ) = cx [U p - 2E k/a ) - co 

= fc" 1 [4(l-a) + (3-2a)Ci]-2C 2 . (B5) 

Finally, we demand that the value of the approxima- 
tion be exact at u p /2. This last condition takes a different 
form dependent on whether u p /2 is greater than or less than 
E k / a as different forms of approximation hold in those two 
regions. Thus 

2 + 2E- 1 {^.) a / k - E- 1 (^) 2/k 
= 2 + dE- 1/a (^ - C 2 E- 2/a {u p /2) 2/k 
for i£ < E k/a , and 
2 a/k E + u p - \u 2 p = (c + Ci^fO ^ for ^ ^ . (B6) 

The four equations (JB4J), (JB5] i) , (Q35] ii) and JB6}, in 
whichever form is relevant, are readily solved for the four 
constants co, ci, Ci and C2. 

Our orbits can now be found in the form, {t/r) k = 
1 + ecos(m0) for h 2 /Ar k ^ E k/2 where m 2 = fc 2 ci, but for 
larger r, we get = 1 + e* cos [m* (<f> + </>*)]. 
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